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9P ' Abstract 

o , 

I We present a continuous time model of maturation and survival, obtained as the limit of a 

compartmental evolution model when the number of compartments tends to infinity. We establish 
in particular an explicit formula for the law of the system output under inhomogeneous killing and 
when the input follows a time-inhomogeneous Poisson process. This approach allows the discussion 
of identifiability issues which are of difficult access for finite compartmental models. The article 
' ends up with an example of application for anticancer drug pharmacokinetics/pharmacodynamics. 
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' Introduction 

> 

■ Compartmental models are widely used in biology and medicine for the modelling of evolution phenom- 
ena. In particular, these models are very usual in the pharmacokinetics/pharmacodynamics analysis 

CsJ I of toxic effects of anticancer drugs on white blood cells. Several mathematical and computational 

■ aspects of such models remain untouched. For instance, at least for catenary chain models, there is 
, no rigorous rule for the choice of the number of compartments, and many identifiability issues are not 
\ elucidated. In the present study, we consider the limit of a discrete family of compartmental models 

when the number of compartments tends to infinity. The obtained limit looks simpler. It provides 
explicit formulas for the mean occupation of the compartment of interest, and allows the discussion of 
important identifiability issues. Roughly speaking, we replace a finite catenary chain of compartments 
with time-inhomogeneous rates by a two compartments time-inhomogeneous model with lag. Much 
^ ■ of the article is devoted to the mathematical derivation and analysis of these models. We hope that 
beyond the mathematical aspects, the main results - illustrated on a simple example at the end of the 
article - may serve the quantitative biologists. 

The mean occupation of the compartment of interest appears as the expectation of an underlying 
stochastic process. This process is a time-inhomogeneous M/M/oo queue, with an explicit biological 
interpretation of its input and output rates. This leads to a nice explicit Binomial-Poisson formula 
for the instantaneous law. In a way, our approach can be seen as a complement and extension of the 
boxcartrain models considered for example in |18l I41j . see also [201 127j and references therein. It is 
also closely linked with "binomial catastrophe models", see for instance [TU] and references therein. 
The novelty is mainly the space-time inhomogeneous killing, the stochastic interpretation in terms of 
Feynman-Kac's formute, and the computation of the output occupation law when the input follows 
a time-inhomogeneous Poisson process. In practice, the formulas that we derive provide a good 
compromise between computer time and numerical errors, without loss in interpretation, as suggested 
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by the example presented at the end of the article. However, pharmacokinetics/pharmacodynamics 
populational aspects are outside the scope of this work. 

Outline of the rest of the article. The rest of the article is organized as follows. 

In Section [H we briefly introduce the time-inhomogeneous deterministic compartmental systems 
with linear rates. In particular, we present two stochastic interpretations in terms of particles. The 
first one is based on interacting time-inhomogeneous M/M/oo queueing systems, whereas the second 
is based on the occupation of independent random walks on the graph of compartments, subject to 
birth and death. These interpretations are at the heart of the results of Section [2] regarding limits of 
catenary chains of compartments. 

In section [21 we present a time-inhomogeneous maturation and survival model, and the related 
counting processes. We show how a finite catenary chain of compartments leads, when the number 
of compartments tends to infinity, to a simple time-inhomogeneous two compartments model with 
lag. The main result is given by Theorem 12.41 which provides simple explicit integral formulas for 
the compartment of interest. These formulas can be seen in turn as a byproduct of Feynman-Kac's 
formulas, and allow the discussion of identifiability issues in Section [2. 3. 4i 

In Section O the results of Sections [T] and [5] are illustrated on a simple example related to anti- 
cancer drug toxicity pharmacokinetics/pharmacodynamics context. Namely, we provide a comparison 
between classical finite catenary chains of compartments in one hand, and our two compartmental 
limit with lag in the other hand. 

1 Finite compartmental systems with time-dependent linear rates 

Consider a system of compartments indexed by the finite set /. Each compartment contains a quantity 
of matter. As we will see in the sequel, the amount of matter can be represented by a discrete 
number, or by a real number, depending on the interpretation chosen. The matter is the subject of 
transfers between compartments. It can also be created (external inflow) or destroyed (outflow) in 
each compartment. 

The reader may find an accessible introduction to simple compartmental systems with several 
examples in the recent book [35] by Matis and Kiffe, see also [26] and references therein. The study 
of the total amount of matter in the system is for instance addressed in |25j , see also |34] . The reader 
may find a study of time-delayed compartmental systems in [37] and [27] and references therein. The 
literature regarding the stochastic interpretations of time-inhomogeneous compartmental systems is 
less rich. Let us recall the essential aspects of finite compartmental systems with linear rates. 

1.1 Deterministic systems and linear ordinary differential equations 

We consider here that the amount of matter is represented by a non-negative real number. Let Qi{t) 
be the amount of matter in compartment i at time t ^ to, where to is the initial time. We denote by 
Q{t) the vector i ^ I \—> Qi{t), and we identify the set R-'^ with M" where n := card(I). Let us consider 
the dynamics of {Q{t))t^to described by the system of linear ordinary differential equations 



with initial condition Q{to). Such a dynamics is described by figured! Here Xi{t) is a creation rate 
for compartment i (inflow), Kj(t) is the destruction rate for compartment j (outflow), and Pij{t) with 
i 7^ J is the transfer rate from compartment i to compartment j. These rates are non- negative. Apart 
from the Aj rates, the rates Kj and pij act proportionally to the content of the compartment that 
gives matter. Note also that a Kj with negative values might mimic a proportional auto-inflow. In 
vector/matrix language, the dynamics of Q{t) is of the form 



(1) 




dtQ{t) = M.{t)Q{t) + \{t), 



(2) 
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where A(i) denote the "vector" (Aj(t); i & I), and where M(t) is the matrix defined by 

Mij{t) :-- 

for any i,j € I. The traditional theory of hnear ordinary differential equations gives 



Q(t) = R(i,to)Q(io) + / B.{t,u)X{u)du 

Jto 

for any t ^ to, where the resolvent R is the solution of the matrix linear differential equation R(fo, ^o) = 
I and 5tR(t, to) = ^o) for any t ^ to- In the literature related to compartmental systems, 

M is sometimes referred as the "transfer matrix". When M does not depend on time, R(f,ii) = 
exp((?; — ii)M). This matrix exponential can be explicitly computed in very special situations where 
for example card(I) is small or where M has a nice structure. Usually, the method is to diagonalize M, 
which leads to an expression of the solution as a linear combination of exponential functions related 
to the spectrum of M. There is no simple closed formula for Q{t) when the rates depend on time. 



Qjit)Pj,i{t) 



Qi{t)pij{t) 




Qi{t)Ki{t) 









Qi{t) 




► 


Pj,iit) 




Ki{t) 



Figure 1: The left hand side diagram shows the flows in a compartmental system. Here, the transfer 
rates and the outflow depends linearly on the content of the compartment that gives matter. In 
contrast, the inflow does not depend on the content of the compartments. It is customary to make 
the Qi{t) and Qj{t) parts of the rates implicit, as show in the right hand side diagram. 



Example 1.1 (Finite catenary chain of compartments). It corresponds to a finite system of n com- 
partments, labelled by I = {1, . . . ,n}, for which Xi = if i ^ 1, and pi j = if j ^ i + 1. It is 
customary to abridge Ai by A and pj,j+i by pi. In such a compartmental system, there is only one 
external inflow with rate X at the left extremity of the chain. Moreover, the interaction of a compart- 
ment is limited to its right neighbor. The topology of this compartmental system is depicted in figure 
[H The matrix M is band diagonal, which makes possible yet tedious the explicit computation of its 
exponential, when the rates do not depend on time. 



1.2 Stochastic interpretation 

Let us consider a finite system of n compartments, labelled by I, and with rates A, k, and p. In the 
stochastic interpretation, each compartment contains particles, and we denote by the number of 
particles in compartment i G / at time t ^ to- The particles are indistinguishable, and can be created, 
destroyed, or can move from one compartment to another, according to a Markovian dynamics of the 
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Figure 2: Finite catenary chain of n compartments with single inflow at the chain head. 



vector Nt := {Nl; i ^ I). For a good choice of Markovian dynamics, the average number E(A't | A't,j = x) 
of particles per compartment is the solution of a system of linear differential equations with initial 
condition x. The reader may find a general presentation of Markov processes and related topics in the 
book [TT] by Ethier & Kurtz. Kelly gave in \6i)\ sec. 4.5 p. 113-117] an accessible introduction to the 
stochastic interpretation of compartmental systems. Here we focus on the time-inhomogeneous case. 

Theorem 1.2. Let I he a finite set. Let {Nt)t^to ^6 a time-inhomogeneous Markov process with state 
space E := N^, and generators (Lt)t^to ■ Assume in addition that for any t ^ t^, the function x £ E 
At{x) := Yliy&E^t,x,yV is affine. For any x £ E, and any t ^ s to, let Q{s,t,x) := E{Nt \ Ng = x). 
Then Q is the solution of the linear differential equation 

dtQ{s,t,x) = At{Q{s,t,x)), 

for any t ^ s ^ to, with initial condition Q{s, s, x) = x. 

Proof. The result is a consequence of Chapman-Kolmogorov equations. Namely, for any t ^ s ^ to and 
any function f : E ^R,let Ps,t{f) be the function E^R defined by Ps,t{f){x) := Hf{Nt) \ = x). 
In particular, P^^sif) = f ■ The Markovianity of the process is captured by the Chapman-Kolmogorov 
equation which writes Ps,uiPu,tif)) = Ps,t{f) ^or any t ^ u ^ s ^ tQ. Recall that acts linearly on 
/ as a matrix. We denote by L^/ the function E ^ R. defined by (Lt/)(x) := J2yeE^t,x,yf{y)- The 
definition of L gives 

This yields, by using the Chapman-Kolmogorov equation, to the "forward equation" 

The desired results follows immediately by considering the ad-hoc / function. Namely, for any i £ L, 
consider the function defined by f{x) = Xi for any x £ E. One has (Ltf){x) = At{x)i. Since At in 
affine and P^,* is Markov, they commute and we get Pt(Ltf){x) = At{Ps,tif){x))i. □ 

Note that the Chapman-Kolmogorov equation which appears in the proof of theorem 11.21 gives also 
the "backward" equation dsPs,t{f){x) = — Ls(Ps^t(/))(a;), which yields the system of linear differential 
equations 

dsQ{s, t,x) = - '^ Ls,x,yQ{s, t, y) 

y&E 

with initial condition Q{s,s,x) = x. Of course, when the process is time- homogeneous, does not 
depend on time anymore and Pg^t = PtoM+t-s- Thus, in that case, Pu(L{f)) = L(Pm(/)) for any 
u ^ to, which makes the forward and backward equations identical up to a sign. 

Let us give now an Lt matrix on E such that At is exactly the affine function of the deterministic 
compartmental system ([2]) constructed from the triplet of rates X{t), K{t) and p{t). We identify the 
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countable space E : — with N'^ where n : — card(/). We denote by 61, • • • , Cn the canonical basis of 
R", embedded in E. For any x,y ^ E with x ^ y, we set 




Xi{t) if y = X + Ci for i G / (birth of a particle) 

XiPij{t) if y = X — Ci + Cj for i 7^ j in / (transfer of a particle) 

XiKi{t) if y = X — ej for i e / (death of a particle) 

otherwise 



The diagonal terms are such that the sum over each row is zero. The transition related to p can 
equivalently expressed as XjPj^i{t) \i y = x — ej + e-i for j 7^ i. Recall that Nt is an n-components 
vector which represents the amount of particles in each of the n compartments at time t. A simple 
computation shows that for any x ^ E and any i £ I, 



We recognize immediately the equation ^ of the deterministic compartmental system. 

1.2.1 Interpretation as interacting M/M/00 queues 

The generator ([3]) suggests an interpretation of the process {Nt)t^to as an n-dimensional M/M/co 
queueing system with interactions. In such a picture, each particle is a client, and each compartment 
is an M/M/00 queue. The quantity is thus the number of customers in the i*'^ queue at time t. For 
each compartment i £ I, and in absence of interaction (p = 0), the arrival rate is Aj, and the service 
rate is Since the queue is M/M/00, each new client gets immediately its own dedicated server, 
which explains the coefficient XjAtj in L. The interaction p allows the clients to move from queue i to 
queue j with rate XiPij. Thus, in presence of interactions, the arrival rate in queue z is + XjPj,i 
and the service rate is XiKi + Xi pij. The clients in the queues are not ordered since these queues 
are of M/M/cxd type. The random vector Nt gives the number of clients in each queue at time t. The 
dynamics of {Nt)t^to can be thus interpreted as n interacting M/M/00 queues, e.g. as n interacting 
birth and death processes on N. These interacting queues can be seen as special Markovian interacting 
particle systems, see [9l[33]. In statistical mechanics, they constitute a "zero-range dynamics on the 
graph of compartments", see [3]. For the reader familiar with queueing systems, these interacting 
queues can be seen as time-inhomogeneous "Jackson networks" [23\ I40j . These queues are "in 
tandem" in the case of a catenary chain of compartments with killing only at the end of the chain. 

In particular, with such an interpretation in mind, a one compartment system is equivalent to a 
single M/M/00 queue. More generally, a catenary chain of n compartments is equivalent to a system 
of n interacting M/M/00 queues, in which only the first one has an external inflow, the others being 
inflowed by the clients that leave the preceding queue due to the interaction. 

1.2.2 Interpretation as independent particles on the graph of compartments 

Following \30\ sec. 4.5 p. 113-117], by removing the indistingishability of the particles assumed in the 
preceding interpretation, Nt is the occupation vector of the compartments at time t for a superposition 
of independent particles moving in /. These particles are subject to birth (A rate) and death (k rate), 
in addition to their motion in the system of compartments viewed as an oriented graph with n vertices. 
This graph is a subset of / x /, and pij = means that the edge does not exist in the graph. In 
other words, p defines an oriented graph, i.e. a topology on the system of compartments. 

Note that when n = 1 (trivial graph), we just have constructed the M/M/00 queue from inde- 
pendent particles with survival rate k, arriving according to an independent Poisson point process of 
intensity A. When n ^ 1 and p = 0, the process consists in n independent M/M/00 queues. 





y&E 
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M/M/oo queue 
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M/M/oo queue 



Random walks 
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Figure 3: The two stochastic interpretations of finite compartmental systems with hnear rates. Inter- 
acting M/M/oo queues versus independent random walks on the graph of compartments subject to 
birth and death. 



In this interpretation, the creation of particles on the n sites is obtained with n independent 
Poisson point processes, one for each external inflow. The independence of the particles moving on 
the graph is a consequence of the linearity in x of the coefficients related to the p rates in hx.y. The 
single particle dynamics involves the p rates for its motions. The k rates are incorporated by the 
addition of an external cemetery state to the system, whereas the A rates are incorporated by using 
an independent Poisson point process and independent copies of the particle. It is known that adding 
an external cemetery state leads to Feynman-Kac type formulas related to a very simple potential, 
see [H] for instance. 

We have thus two stochastic interpretations of the dynamics ([3]) of (Nt), as presented by figure El 
The first one is "vertical" and corresponds to interacting birth and death processes in N (M/M/oo 
queues subject to interactions), whereas the second one is "horizontal" and corresponds to independent 
particles moving on the graph of compartments (random walks on / subject to birth and death). Both 
hold for time-dependent rates. 

2 A continuous stochastic maturation model 

In the sequel, we consider a time-inhomogeneous stochastic maturation model with killing, which 
provides a Poisson-Binomial formula for the occupation law of the maturation system. The probability 
of survival in such a model can be for instance the probability of traversal of a finite catenary chain 
of compartments. We show that this probability has a nice explicit formula when the length of the 
chain tends to infinity. Such a result is in accordance with the stochastic interpretation of finite 
compartmental systems presented in Section [H in terms of time-inhomogeneous M/M/oo queues. 

2.1 Counting mature particles under killing and Poisson production 

Consider a maturation system in which the maturation takes a deterministic duration r > 0. A 
particle which begins its maturation at time T € M achieves its maturation at time T-\-t. We say that 
T + T is the maturation time of the particle. Assume in addition that the particle can die during the 
maturation process, with probability l—p{T), where p(T) G [0, 1] is a deterministic real number which 
may depend on T. In such a maturation with killing scheme, a particle which begins its maturation at 
time T G M can either die with probability 1 — p{T), or survive with probability p{T) and achieve its 
maturation at time T + t. Now, consider independent particles which begin their maturation at times 
To,Ti, . . . These particles survive with probabilities p{Tq),p{Ti), . . ., and achieve their maturation at 
times Mq = Tq + t, Mi = Ti + t, . . . The following Lemma is a consequence of the stability by thinning 
of Poisson point processes, see for instance [31] or [19]. 



6 



Lemma 2.1. Suppose that Tq,Ti, . . . are random and distributed according to a Poisson process with 
intensity A : M — > M+. Assume that this random process is independent from the maturation process 
of the particles. The maturation times Mq,Mi, . . . follow a Poisson point process on M with intensity 
A* : M ^ M+ given by 

X*{t) := X{t-T)p{t-T). 

In other words, the number of mature particles produced during a bounded time interval / C M 
follows a Poisson distribution with mean 

A* {u) du = A(u — r) p{u — r) du. 
I J I 

Moreover, for any disjoint and bounded time intervals I and J, the number Nj and Nj of mature 
particles produced during these time intervals are independent Poisson random variables. The reader 
familiar with point processes may notice that the resulting point process is thus an inhomogeneous 
marked Poisson point process, with Bernoulli marks. 

Let us assume now in addition that independently, each mature particle survives after maturation 
during a random duration with intensity ^ : M+ M+. In other words, given that a mature particle 
is still alive at time s, and if S is its remaining survival duration, then P(S' > t) = exp(— f^fi{w) dw) 
for any t ^ s. The following Theorem expresses the law of the number Nt of mature particles still 
alive at time t € M+. Figure H] gives a schematic view of the system. 

The following Theorem is a direct consequence of Lemma 12.11 but can also be seen alternatively 
as a special case of a result of Keilson &: Servi on time-inhomogeneous M/G/cx3 queuing processes 
(see also [lU EH US] for related problems). Here the resulting counting process {Nt)t^o is a time- 
inhomogeneous M/M/oo queue. The reader may find a general introduction to Queuing Processes in 
many books such as [H E [El EQI ED . 

Theorem 2.2. For any times ^ s ^ t and any m G N, 

CiNt \Ns=m)= B{m, a{s, t)) * V{(3{s, t)) 

where ^ ^ 

a{s , t) := ex.p ^— J fi{w) dw^ and P{s,t) := J X{u — t)p{u — T)a{u,t) du. 

In particular, 'E{Nt \ Ns = m) = ma{s, t) + f3{s, t) for any ^ s ^ t and any m € N. 

Note that there are three sources of randomness here, which are supposed independent: the initial 
time of maturation, the survival during the maturation, and the survival duration after maturation. 



2.2 Modelling the maturation process itself for one particle 

Let us consider a single particle which begins its maturation at time T G M. This section is devoted 
to the modelling of the maturation and survival process itself, by means of a Markovian dynamics on 
the time interval [T, +oo). In particular, it provides an explicit formula of the probability of survival 
during maturation p{T). 



2.2.1 Finite state model 

We begin by considering a maturation model which consists in n + 1 maturation steps, with possible 
killing. More precisely, define the finite set 

:={0,l,...,n}U{c} 

where c {0, 1, . . . ,n} will serve as a cemetery state. The maturation of the particle is modelled by 
a motion from state to state n. State n corresponds to "maturity". The killing rate is modelled by 
a function 

: M X {0, . . . , ?i} ^ M+ 
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Input 



time T 




Output 
time T + r 



Maturity 



Elimination with probability 1 — p{T) Elimination with rate 

Figure 4: Schematic view of a maturation system with killing. Recall that r is the maturation duration. 
The quantity p{T) is the survival probability of a particle which has begun its maturation at time 
T. The quantity (i is the survival rate of a mature particle. The left hand side compartment (LC) 
corresponds to a Poisson deformation box with lag, whereas the right hand side compartment (RC) 
corresponds to an M/M/oo-like counting process. The output of the LC is a point process which 
serves as an input for the RC. 



which is smooth in the first variable (time). This killing rate may thus vary in time and in space. The 
maturation rate is modelled by a constant positive real number p^^\ Consider now the continuous 
time Markov process {x[^^)t>T with state space S^'^\ and generator given for any i 7^ j in 5'^"-' and 
any t ^ T by 







if z 7^ c and j 
if i 7^ c and j 
otherwise 



i + 1 
c 



(5) 



The diagonal terms of the matrix L^") are such that the sum over each row is 0. The dynamics of the 
maturation process {X^^^)t^T can be read directly on the expression of the generator above. Namely, 
at time t and from state i, the process can only move to state i + 1 with rate p^'^\ or die (being killed) 
with rate K,{t, i) and placed in the cemetery state c from which it cannot escape. The maturation steps 
are depicted by figure O In absence of killing, i.e. when k^") = 0, the process {xj;'^^)t^T is a simple 
Poisson process of intensity p^'"^ , stopped when it reaches state n. 



Jn) 



(t,0) 



p 



(n) 



Jn) 



Jn) 



P' 



(n) 



maturity 



flit) 



Jn) 



it,n) 



Figure 5: A catenary chain of compartments with one-way fiow and possible killing, used as a matu- 
ration model maturation (or evolution) rate p^"'\ and killing rate k^'^\ The down arrows lead to the 
cemetery state, which is not represented here. 



In such a maturation with killing model, we are interested in the mean occupation of the state n for 
a Poissonian number of such processes evolving independently. One can define p{T) as the probability 
of hitting state n starting from state 0: 

p{T) = F{3t > T; ^ = n \ X^"^ = 0) and p{t) = /«(")(*, n). 

However, the corresponding maturation time is random, depends on T, and is given by 

T = inf {t > T; X^"^ = n} 
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on the event {x!^^ = 0,3t > T; Xj"^ = 7i}. We will see in the sequel that letting n tends to oo provides 
a very simple explicit formula for p{T) and makes r deterministic, without any loss in interpretation. 

2.2.2 Continuous state limit model 

It is possible to express the law of (X("))t^T by mean of a Feynman-Kac formula. Rather than 
using this approach, we prefer to consider now the limit process obtained by letting the number of 
compartments n tends to oo. We will use a scaling which corresponds to the convergence 

1 n 



; ) ' ' ' 5 

n n n 



.[0,1]. 

Namely, assume that p^"^ = np for some positive constant p. Assume that 



lim sup 

"^°°iG{0,...,n}; 



- K{t,i/n) 



for some bounded piecewise continuous function 

k:Mx [0, 1] 

Theorem 2.3. Let x € [0,1]. On the event {x!p'^ = \nx\}, and with the convention c/n = c, the 

resettled process {n~^ x[^^)t^T converges in Ittw, when n goes to oo, toward the Markov process {Xt)t^T 
with state space S := [0, 1] U {c} defined for any T ^ s ^ t by 



jO{Xt I Xs 



p{s, t, x) (5^(f_^) + (1 - pis, t, x)) 6c ifx^c 
5 c if X = c 



where x{u) := min(l, x + p{u — T)), p{u, v, x) := exp (— /J k{uj, x{w — u)) dvj) , and p{u, v, c) := for 
any x / c. 

Proof. Consider the natural inclusion 7r„ : S^^^ S defined by 7r„(i) = i/n for any i G {0, . . . , n} and 
^n(c) = c. By this way, any function / : 5 ^ M induces naturally a function /(vr„) from 5^"^ to M, 
and one has 

L-(/W)(.)^n.(/(i±i)-/(i)).„(.i)(/(c,-/(i)). 



By following step by step the classical argument based on Taylor's formulas presented in |1H p. 1-5], 
one can show that the Markov process (n~^ x'f^^)t^T converges in law, when n goes to oo, toward the 
Markov process {Xt)t^T with state space S := [0, 1] U {c} and infinitesimal generators 



I if oc — c 

defined for continuous functions / : 5 ^ M which are smooth on [0, 1] and vanish at the boundaries. 
The scaling method that we used is related to the fluid limits of Queuing Theory (see for instance 
in |40t ch. 9]) and the hydrodynamic limits of interacting particles systems in Statistical Mechanics 
(see for instance (32]). The obtained generator above is the addition of a deterministic constant 
drift term pf'{x) together with a random space-time inhomogeneous killing term K{t,x)(f(c) — f{x)). 
A Feynman-Kac formula or a direct check gives the desired expression of the law. Actually, the 
addition of a cemetery state corresponding to a killing term to a Markov process always gives rise to 
a Feynman-Kac formulae, see In some ways, our case ^ is the simplest one: killing term + 

deterministic generator (constant drift). 

□ 
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The interpretation in terms of the position Xt of the particle is quite clear: the particle is moving 
from left to right on the interval of positions [0, 1] at constant speed p. The particle stops its motion 
when (and if) it reaches the right extremity of [0, 1]. At any time t G M+ and any position x E [0, 1], 
it can be killed (and thus placed in the cemetery state c) with rate K{t,x). A particle is mature when 
it reaches the right extremity 1 of [0, 1]. The positions [0, 1) correspond to the steps of maturation. 
For any t € M and any x G [0, 1], let us define g{t,x) := K(t, x)I[o^i) (x) and n{t) := nit, 1) in such a 
way that 

K{t,x) = 5r(t,x)I[o,l)(x) +/i(t)I|i}(x). (7) 

Function g gathers the killing rate during maturation, whereas function p captures the killing rate 
after maturation. In the settings of Section 12.11 the process {Xt)t^T corresponds to a maturation 
model for which 

fi{t) = K{t,l) and p{T) = p{T,T + T,0) = exp ^— J g{T + w, pw) dw^ 



pw) dw ] , (8) 



where the "maturation time" r is deterministic, independent of t, and given by 

.:=i. (9) 

When the particle is not killed, i.e. on the event {Xt = 0,3s > 0;Xs = 1}, r is exactly the 
deterministic duration taken by the particle to reach state 1 from state (at constant speed p). 

2.3 The resulting counting process 

Consider the maturation system with killing modelled by the {Xt) process introduced above, with 
maturation speed p and killing rate k : M x [0, 1] M4. as in ([Tj). The maturation takes a deterministic 
time r given by ([9]). A particle which begins its maturation at time T can either die with probability 
1 — p{T) where p{T) is as in ([H]), or survive with probability p{T) and achieve its maturation at time 
T + T. Now, consider independent particles which begin their maturation at times To,Ti, . . . These 
particles survive with probabilities p(ro),p(Ti), . . ., and achieve their maturation at times Mq = 
To + r. Ml = Ti + T, . . . provided that the corresponding p(Tj) are non null. The survival durations 
after maturation of the mature particles are i.i.d. random variables with common rate p : M_|_ — s- M_|_ 
where p{t) := K{t, 1). Suppose that the initial times Tq,Ti, . . . are random and distributed according 
to an independent Poisson point process on M with intensity A : M ^ M+. For any t ^ 0, let Nt be 
the number of mature particles still alive at time t. Theorem 12.21 together with ([8]) and ([9]), yields the 
following result. 

Theorem 2.4. For any ^ s ^ t and any m € N, C{Nt \ Ng = m) = B{m, a{s, t)) * V{(3{s, t)), where 



a{s, t) := exp ^— J 



and 



p{w) dw 



P{s,t) := J X{u)ex.py — J g{u + w, pw) dw j a{u + T,t) du. 

In particular, for any ^ s ^ t and any m G N, 

B{Nt\Ns = m) = ma{s,t)+l3{s,t). (10) 

The process {Nt)t^o is a time inhomogeneous M/M/00 queue, which is a particular birth and death 
process, see |401 [21 HI]. In an M/M/00 queue, each client is immediately served by an independent 
dedicated server. When g = and p is constant, one has p ^ 1, and {Nt)t^o is in that case an M/M/00 
queue with input intensity A and constant output intensity p. When in addition A is constant. 
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the symmetric invariant measure of this queue is the Poisson law V{X/fi). In a way, the M/M/oo 
queue with constant intensities is for the Poisson process what the Ornstein-Uhlenbeck process is for 
Brownian motion, see [40 ^ Theorem 6.14]. 

The role played by the M/M/oo queueing processes in our model is due to the independence of the 
particles in the definition of {Nt)t^o- One can alternatively consider a non independent killing after 
maturation, which can lead for example to an M/M/1 type queueing process. Unfortunately, the law 
at fixed time of such a process, even for the time-homogeneous case, is far more complex than the 
simple formula obtained in the M/M/oo case, see |42] and |40j . 

Remark 2.5 (Negative values of k and input rate amplification heuristics). The expression of (3 
above still makes sense even when k takes negative values on [0,1). In that case, the quantity p{u) = 
exp {^— g(u + w, pw) dw) may exceed 1, and thus cannot be interpreted as a "survival probability". 
Such negative values can be seen in a way as a sort of "amplification" of the input rate instead of 
a killing, and can be incorporated into A. Beware that it does not correspond to an immigration of 
particles in the counting process. It appears as a distortion of the input rate of the counting process. 



Poisson input 



Maturity 






p 


1 


p 


p 


n 











K(t,0) 



Hi{t, 1) 



K[t, n) 



elimination (cemetery state c) 



Figure 6: Catenary chain of compartments with Poisson input, one-way flow (p), and space-time- 
inhomogeneous killing (k). Its limit when n tends to infinity is given by figure [H 



2.3.1 Case without killing after maturation 

Assume that = (no killing after maturation). In that case, a{s,t) = for any ^ s ^ t, and the 
formula for (]{s, t) boils down to 

f5{s,t) = J A(u)exp^— y g{u + w, pw) dvj^ du. 

In that case, t S f3{s, t) is non decreasing since the main integrand does not depend on t. This 

is not surprising since the particles are never killed after maturation. Hence, on {A^^ = 0}, the process 
t € M+ I— > Nt is non decreasing, and in particular its expected value /3(s, t) is non decreasing too. 
When A is constant and g = 0, we recover a simple Poisson process of intensity A. 
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2.3.2 Case without killing during maturation 

Assume that g = (no killing during maturation). In that case, the formula for (3{s,t) for ^ s ^ f 
boils down to 

rt—r 

P{s,t) = / X{u)a{u + T,t) du. 

J S — T 

When in addition both A and are constant, we recover the standard time-homogeneous M/M/oo 
queue with input rate A and service rate /i, for which 

a{s, t) = e-^*-^)/^ and ^{s, t) = [l - e~^'"^^)j^- (H) 

The Poisson measure 'P(A//i) is a stationary distribution of {Nt)t^Q in that case. 



2.3.3 Partial killing during maturation and constant killing after maturation 

Let us consider now the particular case where A is constant, and k is of the form 

K{t,x) = 9(t)l[o,5)(x) +^I{i}(x), (12) 

where 6 G [0, 1], where /i G M-|_, and where g : M-)_ M_|_ is a smooth function. It corresponds to a time 
dependent killing before the maturation stage 5, and to a constant killing after maturation. No killing 
is made between maturation stage 5 and full maturation. The formula for (3(s, t) when ^ s ^ t boils 
down to ^ ^ 

l3{s,t) = X J exp ^— fj,{t — T — u) — J g{u + w) dvj^ du. (13) 

When g = 0, this formula reduces to the classical M/M/oo average queue length (jlip . Assume instead 
that function g in ()12p only vanishes at infinity. Then the two formulas (1131) and (jlll) for /3 are 
equivalent when t goes to +00. In particular, the Poisson measure V{X/ ^) is a stationary distribution 
of {Nt)t^o in that case. One can see on figure [7] that this Poisson equilibrium is quickly reached. It is 
possible to quantify the speed of convergence in total variation norm or in entropy, see [U HQ] . 



2.3.4 Identifiability and invariance by some transformations 

The dynamics of N := {Nt)t^o is fully described by the quadruple {T,iJ,,g,X), where r := 1/p. It is 
quite natural to ask about the injectivity of the map 

iT,fi,g,X)^£iN\No). 

According to Theorem 12.41 the law £(A^| A'^o) is completely prescribed by the triple (r, /i, A^), where 
Ag : M ^ M+ is defined by 

Xg{t) := X{t)p{t) = X{t) exp j g{t + w, pw) dw^ . 

Thus, the couple {X,g) cannot be identified, since it acts on the dynamics via the compound A^,. 
The action of g can be compensated by A and vice versa. Namely, suppose that g can be written 
g = gi + g2, where gi and g2 are non-negative functions. Then, the two models corresponding 
respectively to (r, /i,g,A) and {t, p., g2, Xg^) are indistinguishable. The extreme case corresponds to 
{t, p,0, Xg), for which the entire killing function g is merged into the input rate function A. Let us 
analyze some special cases. For any ^ > 1, let us consider the transformation which replaces {X,g) by 
{X'^,g'^) defined by 

A^ := ex and g^ ■= g + p log{9). 

Function p and parameter p are left unchanged, and one can check that A^g = Xg. Therefore, the 
dynamics is invariant by this transformation: the models corresponding to (r, p, g, A) and (r, p, g^ , A^) 
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Figure 7: The plots represent the average number K{Nt \ Nq = 0) of Theorem 12.41 the killing 
function g. This example corresponds to the case (fT^ with A(t) = iR^(t), a = ^, // = 1, /) = 1, and 
g{t — 10) = (e"*/^'' — e~*/^)I]R^ (t). One can observe on the plot of the average number of particles 
three main time phases. The first phase corresponds to an ascendancy to a Poisson equilibrium before 
the action of the drug via g. The second phase corresponds to a decrease due to the action of the 
killing via function g (delayed by the time lag is r = = 1)- In the third phase, the killing action 
decreases and the Poisson equilibrium is reached again. 



are indistinguishable, despite their distinct physiological meanings. A multiplicative perturbation 
of the input intensity A corresponds to an additive perturbation of the killing function g during 
maturation. Hence, one can decide to normalize the parametrization by taking for example A = 1. 
Namely, (r, fi, g, A) is indistinguishable from (r, /i, G, 1) where 

G{v, y) = g{v, y) - - log X ( v - -] . 



3 A pharmacokinetics/pharmacodynamics example in cancerology 

Catenary chains of compartments, as depicted in figured are widely used in the literature by kineticists 
for the modelling of anticancer drug toxicity, see for instance [6l[7tll4 t ll5 t ll6 t [T7] and references therein. 
Most anticancer drugs have toxic effects on white blood cells (myelosuppression) . Neutrophils are 
particular white blood cells used as toxicity markers. The catenary chain of compartments models 
the maturation process of neutrophils in the bone marrow. Each maturation stage corresponds to a 
specific position in the chain. The killing rate during maturation corresponds to the drug toxicity on 
neutrophils. The last compartment of the chain is the only observed (blood), and the others correspond 
to hidden positions (bone marrow). In practice, the input and transit rates in the catenary chain of 
compartments are unknown parameters. The goal of the kineticist is thus to control the content of the 
last compartment by his action on the anticancer drug dosage regimen. Consider for simplicity that 
both the production rate A, the transit rate p, and the elimination rate in blood /i are constant. The 
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killing rate k depends typically on the time profile q : ^ M of the anticancer drug concentration 
in blood (i.e. kinetics). The first step of the Pharmacokinetics/Pharmacodynamics study is to find, 
for a given kinetics, the best values of the parameters in view of the data (in blood). This first step 
can be addressed in several ways. 

Approach (I) . This traditional approach involves a deterministic finite catenary chain of n compart- 
ments, as presented in example 11.11 page [3l It is customary to model the toxic effect of the anticancer 
drug on neutrophils via a killing rate K(t, i) = ^q{t) for i ^ hq < n, K{t, n) = fi, and elsewhere. The 
vector Q{t) := (Qi{t), . . . , Quit)), which represents the number of neutrophils in each compartments 
of the chain, solves a system of n ordinary differential equations (ODE). The available data concerns 
only the last compartment (blood), labelled n. The function t i— > Qn{t) depends on the modelling 
parameters A, p, fi, 7, no. The main drawback of this approach is that the number of compartments 
n is unknown, and that the system of equation is heavy. 

Approach (II). In this approach, we circumvent the problem of the choice of n by considering the 
continuous limit in n, as presented in Section 11.21 This limit procedure leads to killing rates k, as in 
equation ([7]) . More precisely, we mimic the cutoff used in the approach (I) by considering equation (fT2|) 
with g{t, x) = 7g(t)I[Q 5] (x) where 7 > and where 5 S [0, 1) are constants. The continuous limit in n 
is constructed at the level of the stochastic interpretation. As a consequence, the relevant quantity is 
the mean of the counting process of Theorem 12.41 given by equation (|1U|) . The initial value m in (|lUp 
must be averaged with respect to the equilibrium without drug, i.e. the Poisson measure V{\/ p). The 
quantity Qn{t) of the approach (I) is thus replaced now by the following explicit alternative quantity: 

goo(t) :=a(0,t)+/3(0,t)-, (14) 

where the functions a and (5 are as in Theorem 12.41 

Remark 3.1. Actually, one can use a Hill-type transform (see for instance ]36^ ) to model the depen- 
dency of the killing rate with respect to the drug kinetics. However, for low drug concentration values, 
it is customary to simplify the approach by considering a linear dependency with a possible cutoff in 
space, as in the approaches (I) and (U) above. We emphasize the fact that even if the instantaneous 
killing rate depends linearly on the drug kinetics, the global killing effect is nonlinear with respect to 
the kinetics. 

We compare below, on a simple example, the two approaches defined above. The approach (I) 
leads to the resolution of a system of ODE, whereas the approach (II) leads to the computation of 
explicit integral transforms given by (I14p . 

In order to mimic practical situations, we used a synthetic dataset obtained by perturbing a dataset 
extracted from [6] (we did not have the permission to use directly the dataset of [6]). In this study, 
an anticancer drug was administrated by 30-min intravenous infusions to patients, on five consecutive 
days, see figure [HI For a typical patient, the drug concentration at time t was proportional to : 

m = E (1 - ^-"^*-'''0^24.+[o,i/2](t) + (1 - e-"/^)e-^(*-2^'^-V2)i,,,^(^/,^^,(t), 

with a = 1.86 and (3 = 0.51. For simplicity, we restrict our analysis in the sequel to a single patient. 
We first analyzed this dataset using approach (I) and approach (II), as described above. Let us denote 
by Yj the observed neutrophils counts in blood at time tj. 

Approach (I). The vector valued function t {Qi{t),... ,Qn{t)) solves the following system of 
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ordinary differential equations: 



' dtQiit) = \-pQi{t)-^q{t)Qi{t) 
dtQ2{t) = pQi{t) - pQ2{t) - -fq{t)Q2{t) 

dtQnoit) = pQno-l{t) - PQnoit) - iq{t)Qno{t) 
, dtQn{t) = pQn-l{t) - P'Qnit) 

At time t = 0, the drug concentration q is null, and the system is at equilibrium. Consequently, 
dt=oQi = for i € {l,...,n}. The initial condition is thus Qi(0) = ••• = Q„_i(0) = A/p, and 
Qn(0) = A//i. These systems were numerically integrated using the Fortran package ODEpack, 
see for instance [Ml EEl ES] • Models with n= 5, 10, 30 and 100 compartments were used. For each 
n-compartments model, the parameters (A, p, p, 7, no) were estimated with ordinary least-squares : 

arg inf y^{Y, - Qn{t,)f . 

3 

Since ng is an integer, all possible values of no, from 1 to n — 1, were screened and the value of no 
giving the minimum residual sum of squares was then selected. 

Approach (II). The explicit expression of (I14p provides the following formula for the average number 
of neutrophils in blood at time t: 

Q^(t) = -e"^* + Ae''^*"^) y exp^pu-jj q{u + w)dw^ du, 

where r = 1/p. Since q has an easily computable primitive, the inner integral in the expression above is 
explicit. The outer integral was numerically evaluated using the Clenshaw-Curtis quadrature method. 
The parameters (A, p, p, 7, 6) were estimated with ordinary least-squares. Note that the continuous 
parameter 5 € [0, 1) plays here the role of the discrete parameter no < n of the approach (I). 

The curves obtained with the different models are represented in figure [HI On the whole, the 
different models approximately give the same curves. Despite its appearance, t 1-^ Qooit) is smooth. 
The curves provided by the approach (I) with 100 compartments in one hand and the approach (II) 
in the other hand are nearly the same. The obtained residuals sum of squares are given in table [TJ 
The curves are given by figure [8l The best fit was obtained with the 10-compartments model. Models 
with more compartments, including the 00 compartments model of approach (II) give approximately 
the same quality of fit. The number of parameter to estimate is the same for both approach (I) and 
(II). However, the estimation of no in approach (I) leads to a large amount of extra ODE integration. 
The amounts of time needed for the ODE integration are at least twice as long for models with more 
than 10 compartments than for the integral evaluation of approach (II). Keeping in mind that such 
studies usually include several tens of patients, the approach (II) that we proposed offers a reasonable 
alternative to compartmental models. 





5 CP 


10 CP 


30 CP 


100 CP 


Continuous 


Computation time 


0.39 


0.99 


11.803 


235.12 


0.51 


Sum of Squares 


0.99 


0.25 


0.32 


0.32 


0.33 



Table 1: The duration of the ODE/integral evaluation. For compartmental models, these durations 
have been obtained as (n — 1) x duration for a single evaluation of the ODE to take into account the 
screening needed to estimate no. The obtained sum-of-squares are given in the second row. 
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Figure 8: The curve with five peaks is the drug concentration curve t i— > q{t) (kinetics). The neutrophils 
counts observed in a patient are represented by squares. The mean neutrophils counts predicted with 
models with 5, 10, 30, 100, oo compartments are represented with different lines. On the whole, models 
with at least ten compartments give a reasonable description of observed points. We emphasize that in 
general, the adequate number of compartments is not known a priori, and can vary from one patient 
to other. The estimation of an adequate number of compartments from the data has a cost. The 
continuous model avoids this problem. 
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